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FOREWORD 

This  report  documents  a  computer  program  which  solves  an  improperly  posed 
heat  conduction  problem  arising  in  the  measurement  of  heat  transfer  data  in  wind 
tunnel  experiments.  The  mathematical  model  which  is  treated  is  one  dimensional 
and  assumes  constant  material  properties.  The  computational  approach  is  based  on 
an  adaptation  of  a  method  used  by  Alifanov  to  treat  a  similar  problem.  The  Green's 
function  for  a  well  posed  heat  conduction  problem  is  used  to  reduce  the  problem  to 
a  Volterra  integral  equation  of  the  first  kind.  Tikhonov  regularization  techniques 
are  then  applied  to  compute  stable  approximate  solutions.  Numerical  results  for 
several  cases  typical  of  a  class  of  wind  tunnel  experiments  are  presented. 


IRA  M.  BLATSTEIN 
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1.  INTRODUCTION 

One  of  the  aerodynamic  properties  to  be  measured  during  a  wind  tunnel  experiment 
is  the  heat  transfer  rate  to  the  surface  of  a  thin  walled  wind  tunnel  model. 
Traditionally  this  is  done  by  assuming  a  functional  form  for  the  temperature 
profiles  and  using  measured  data  to  fit  parameters.  Such  approximations  introduce 
significant  error  in  short  duration  wind  tunnel  tests  which  feature  rapid  model 
pitching.  Unfortunately,  improving  upon  this  simple  approximation  is  not  straight¬ 
forward  because  the  mathematical  problem  describing  the  relationship  between  the 
measured  data  and  the  desired  information  is  improperly  posed  and  is  consequently 
unstable  with  respect  to  small  measurement  errors. 

Various  treatments  of  ill-posed  heat  conduction  problems  have  appeared  in  the 
literature.  An  excellent  survey  of  engineering  methods  can  be  found  in  Beck.^ 
Somewhat  more  mathematical  approaches  can  be  found,  for  example,  in  Cannon  and 
Douglas ^  and  in  Ginsberg.-^ 

The  approach  used  in  this  report  is  a  modification  of  a  method  used  by  Alifanov4 
to  treat  a  similar  problem.  Essentially,  the  method  recasts  the  problem  as  an 
integral  equation  of  the  first  kind  which  is  ill-posed.  A  Tikhonov  regularization 
procedure^  is  then  used  to  compute  stable  approximate  solutions  to  the  integral 
equation. 

In  the  next  section,  the  mathematical  model  is  presented  and  the  integral  form 
of  the  solution  is  derived  from  the  Green's  function  for  the  heat  equation.  The 
application  of  the  Tikhonov  regularization  procedure  and  the  subsequent 
discretizations  are  discussed  in  Section  3.  The  last  section  presents  some  numeri¬ 
cal  results  characteristic  of  wind  tunnel  experiments.  A  listing  of  a  program 
implementing  the  method  is  given  in  Appendix  A. 
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2.  MATHEMATICAL  FORMULATION 

To  determine  heat  transfer  rates  on  a  body  in  a  wind  tunnel  experiment  it  is 
necessary  to  measure  heat  flux  on  the  surface  of  the  model.  However,  it  is  often 
inconvenient  or  excessively  costly  to  place  thermocouples  directly  on  the  surface 
of  the  model  which  do  not  effect  flow  characteristics.  Instead  it  is  common 
practice  to  place  thermocouples  on  the  inside  surfaces  of  thin  walled  models  which 
results  in  a  direct  measurements  of  the  inner  surface  temperature.  If  we  assume 
constant  material  properties  and  that  heat  conduction  along  the  body  surface  can 
be  Ignored,  the  mathematical  problem  of  determining  heat  transfer  rates  on  the 
front  surface  reduces  to  determining 

<  ux,  a,t*) 

given  that  the  temperature  U  (x',  t')  satisfies 

pc  Ut,  “  k  U^x,  0  <  t'  £  T,  0  <  x’  <  l 


U(0,t*)  -  f(t')  0  <  t'  <_  T 

U(x\0)  -  0  0  <  x’  <  l 

<  Uj(0,t)  -  0  0  <  t'  £  T 

Here  t'  denotes  time  and  x*  corresponds  to  the  distance  from  the  backface  with 
x'  *  0,£  representing  the  inner  and  outer  surfaces  respectively.  Also,  pc  is  the 

A 

volumetric  heat  capacity  and  <  is  the  thermal  diffusivity.  The  function  f(t’) 
respresents  the  measured  temperature  data. 

The  two  additional  boundary  conditions  U(x,0)  ■  0,  k  11^(0, t)  ■  0  are 
reasonable  approximations  of  typical  wind  tunnel  conditions.  The  first  assumes 
that  the  initial  temperature  is  constant  (normalized  to  zero)  while  the  second 
presumes  an  adiabatic  inner  surface.  This  second  assumption  is  reasonable  because 
of  the  extremely  low  pressure  inside  the  model  during  a  test  and  the  short  duration 
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of  the  tests  under  consideration.  Note,  however,  that  the  more  general  case  where 
u(x,0)  and  <  11^(0, t)  are  known  functions  can  be  reduced  to  the  above  form  by  solved 
a  well  posed  auxiliary  heat  conduction  problem. 

The  nondimens ional  form  of  the  problem  obtained  by  setting  t  =*  t'/T,  a  =  ocp/tc 


and  x  *  x'/l  is  to  find 
a  Ux(l,t) 

given  that  U  satisfies 


U  -  a  U 
t  xx 


U(x,0)  =*  0 


U(0, t)  -  f(t) 


ux(0,t)  -  0 


0  <  t  <  1 


0<t<l,  0  <  x  <  1 


0  <  X  <  1 


0  <  t  <  1 


0  <  t  <  1 


In  applying  the  Tikhonov  regularization  techniques  to  the  problem  (1)  it  is 
helpful  to  recast  the  problem  as  an  integral  equation.  We  suppose  that  we  are 
given  a  Ux(l,t)  -  g(t)  and  consider  the  problem 

U_  =  a  U 
t  xx 

U(x,0)  -  0 

a 

u(0,t)  *  0 

a  Ux(l,t)  -  g(t) 

The  solution  to  this  problem  can  be  written  (cf.  Ref.  6)  in  the  form 


U(x,t)  -  fo  0(x , t  — t)  g(T)  dT 


where 


0(x,t) 


,  l"-(x+2n)2  *1 

la  X  SXP  "^t - 

■  ✓rat  L  -i 
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Using  this  integral  representation  the  problem  (1)  reduces  to  finding  g(t) 
such  chat 


t 

f(t)  -  /  0(0, t-T)  g(T)dT  =  Kg  (2) 

o 

Unfortunately  the  integral  equation  (2)  for  g;  reflecting  the  ill-posed  nature  of 
the  original  problem,  is  of  first  kind;  consequently,  small  perturbations  in  f 

A 

cause  large  changes  in  the  corresponding  g(t).  In  fact,  since  f  contains  measure¬ 
ment  errors,  it  will  not  be  smooth  and,  in  general,  there  will  be  no  g  e 
(square  integrable)  with  satisfies  the  operator  equation 

Kg  “  f . 

Nevertheless,  special  techniques  can  be  used  to  compute  stable  approximate 
solutions  to  be  integral  equation. 
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3.  REGULARIZATION  AND  NUMERICAL  APPROXIMATION 

In  Che  previous  section  we  used  an  integral  representation  of  the  solution  of 
(1)  to  reduce  the  problem  to  a  first  kind  integral  equation  of  the  form  (2).  We 
also  noted  that  measurement  error  leads  to  possible  nonexistence  and  instability. 

To  overcome  the  nonexistence  problem  we  note  that  for  the  correct  g  there  is  an  f 
for  which 

Kg  *  f. 

We  now  assume  that  we  have  a  bound  e  on  the  error  introduced  by  measurement  of  the 
form 

fl  (f(t)  -  f(t))2  dt  <  e2  (3) 

o 

A 

Based  on  the  consideration  that  the  true  temperature  profile  at  x  *  0  is  near  f  we 
seek  a  smooth  g  (g  e  H^,  i.e.  g  and  g'  are  square  integrable)  which  minimizes 

| |Kg  -  f | |l2  (4) 

However,  the  least  squares  form  (4)  still  suffers  instability  properties. 

To  stabilize  the  least  squares  problem  (4)  we  use  regularization  as  described  in 
Tikhonov  and  Arsenin.5  The  idea  behind  regularization  is  to  add  a  penalty  to  the 
least  squares  problem  based  on  a  norm  of  g.  In  this  case  we  use  H*  -  regularization 
of  the  form 

min  ||Kg  -  f | |^2  +8  ||g|||l  (5) 

where  8  is  known  as  the  regularization  parameter.  According  to  theory  developed 
in  [5],  if  8  *  8(e)  is  choosen  so  that  the  minimizer  gg  of  (5)  satisfies 

C2  e  <  |  | Kgg  -  f  |  Il2  <  Cj  e  (6) 

for  two  constants  and  C2  then 

gg  -*•  a  ux(1’c>  as  z  "*■ 

Thus,  the  computation  of  stable  approximate  solutions  to  the  integral  equation  (2) 
reduces  to  the  minimization  of  (5),  suitably  adjusting  8  to  satisfy  (6). 

To  discretize  the  minimization  (5)  we  will  restrict  our  consideration  to  a 
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finite  dimensional  set  of  functions  S.  We  will  then  define  a  method  of  associating 

a 

a  function  in  S  with  Kg  for  g  e  S.  Since  f  is  only  given  by  discrete  measured 
-  slues,  the  most  natural  choice  for  S  is  the  piecewise  linear  functions  on  a  given 
grid.  More  precisely,  we  subdivide  the  time  interval  [0,1]  into  N  subintervals  of 
length  At  -  1/N  and  let  S  be  the  set  of  functions  which  are  continuous  and  linear 
on  each  subinterval.  (The  program  assumes  that  N  corresponds  to  the  measured  data 

A  A 

points  for  f  and  that  f  is  the  piecewise  linear  interpolant  of  these  measured  values) . 

The  minimization  (5)  has  now  become  finite  dimensional;  but,  we  have  yet  to 
specify  how  to  associate  a  function  in  S  to  Kg.  One  possible  approach  to  this 
problem  would  be  to  approximate  the  integral  in  (2)  by  a  quadrature  rule  and  use 
this  to  define  nodal  values  for  Kg.  Unfortunately  the  expression  for  0  as  an 
infinite  sum  makes  this  a  cumbersome  task. 

A  better  approach  is  to  recall  that  Kg  is  u(0,t)  for  the  solution  to  the  heat 
conduction  problem  (1').  Using  this  property  and  the  linearity  of  K  allows  one  to 
use  a  numerical  solution  procedure  for  (1*)  to  find  approximate  nodal  values  for  Kg. 
More  precisely  we  define  gg  e  S  to  be  the  piecewise  linear  function  such  that 

g5  (jAt)  *  1  if  j  =  1 
0  if  j  i  1 

Introducing  a  spatial  finite  difference  grid  and  using  the  Crank-Nicholson 
integration  scheme  (cf.  Ref.  7)  we  can  solve  (!')  numerically  with  g  =  to  obtain  a 
vector  of  values  corresponding  to  the  values  of  u  (0,iAt).  The  linear  interpolant 
f^  of  the  nodal  values  f ^  then  approximates  Kgg.  For  most  applications  the  error 
introduced  by  the  numerical  integration  procedure  is  much  smaller  then  £  and  can 
be  ignored. 

We  now  denote  the  components  of  f^  by  f ^  and  define  the  matrix 
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Using  the  linearity  of  K,  superposition  and  denoting  by  g  the  vector  of  nodal  values 
of  g,  we  see  that  Ag  defines  a  vector  which  approximates  nodal  values  of  Kg. 

To  complete  the  discretization  we  note  that  for  g  e  S  we  can  evaluate  the  hA  - 
norm  to  obtain 

1 lal I2  5  /‘  |g(t)|2  +  |g'(r)|2  dt 

0  ' 


The  minimization  problem  (5)  now  becomes 

min  (A^  -  ?)T  (Ag  -  1)  +  3  B  g.  (7) 

-v  T  -v 

(Note  that  (Ag  -  f)  (Ag  -  f)  is  not  the  exact  integral  of  the  associated  piecewise 
linear  function;  but,  is  only  the  sum  of  the  squares  of  the  residuals.  This 
reflects  the  type  of  estimate  available  for  e  whereas  a  higher  order  integration 
would  only  integrate  the  noise  more  precisely.  The  minimizer  of  £g  of  (7) 
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icisfies  che  Euler  equations;  i.e., 

(AT  A  +  SB)  g3  -  AT  f  (8) 

Lnce  A^  A  +  3b  is  symmetric  and  positive  definite,  (8)  is  uniquely  solvable  and 
j  can  be  found  by  Gaussian  elimination  (cf.  Ref.  8). 

The  value  of  3  remains  to  be  specified.  The  theoretical  basis  for  this 
slection  is  che  discrepancy  principle  [5]  of  setting  3  according  to  (6). 

\  practice,  since  the  residual 

r(3)  -  j^(Agg  -  f)'  (Agg  -  f)  J  ^ 

3  a  montone  (increasing)  function  of  3,  3  can  be  set  by  trial  and  error, 
icking  3  so  that  r(S)  was  very  close  to  t  was  found  to  yield  satisfactory  results, 
umerical  experiments  indicated  that  once  an  appropriate  3  was  selected,  it  need 
ot  be  changed  unless  the  character  of  the  noise  changes. 
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NUMERICAL  RESULTS 

To  test  the  algorithm  three  cases  typical  of  a  class  of  wind  tunnel  experiments 
were  treated.  In  each  case  the  length  of  the  test  was  assumed  to  be  1.25  seconds 
and  the  material  was  assumed  to  be  .054  inch  thick  steel  with  thermal  conductivity 
k  »  0.00289  Btu/ft.  sec.  °F  and  volumetric  heat  capacity  pc  *  58.98  Btu/ft^  °F.  To 
obtain  data,  three  heat  flux  distributions  were  specified  and  the  heat  conduction 
equations  were  used  to  compute  numerically  the  temperatures  on  the  back  face 
(x  =  0) .  It  should  be  noted  that  the  assumption  of  linearity,  one  dimensionality 
and  insulated  boundary  at  x  *  0  were  assumed  in  obtaining  the  data;  consequently, 
these  numerical  experiments  do  not  test  the  validity  of  these  assumptions.  The 
specified  heat  fluxes  and  associated  temperature  profiles  at  x  =  0  are  shown  in 
Figures  la, lb  and  lc. 

We  now  have  three  "true"  backface  temperature  profiles  for  which  we  know  the 

associated  heat  fluxes.  The  numerical  test  will  now  be  to  add  noise  to  the  backface 

temperature  distributions  and  attempt  to  recover  the  associated  heat  fluxes.  Two 
noise  levels  were  used.  They  are  shown  in  Figures  2a  and  2b.  The  first  is  typical 
of  wind  tunnel  thermocouple  measurement  noise  and  the  second  prescribes  a  noise 
level  much  higher  than  expected  in  practice. 

For  the  low  noise  level  case,  it  was  found  that  for  8  *  10“^  the  computed 
residual  rms  noise  level  agreed  well  with  the  true  rms  noise  level.  The  numerical 
results  for  the  three  cases  are  shown  in  Figures  3a,  3b  and  3c.  The  true  solution 

is  represented  by  the  solid  line  and  the  computed  solutions  is  shown  by  the  dotted 

line . 

Note  that  in  each  case  there  is  poor  agreement  between  true  and  computed 
solutions  neat:  the  end  of  the  interval.  The  cause  of  this  problem  is  that  the  heat 
flux  near  the  final  time  has  very  little  effect  on  temperatures  at  x  *  0.  For  this 
reason  the  second  term  in  (5)  dominates  the  minimization.  The  dominant  term  than 
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UCUUUiCa  /A  -  ,7 

/  |g *  (t>  r  dT 

o 

so  that  g'  (t)-»  0  as  c  -*•  T. 

For  the  high  noise  case  the  appropriate  value  of  the  regularization  parameter 
was  found  to  be  8  ■  .001.  In  Figures  4a,  4b  and  4c  we  show,  for  each  case,  the 
noisy  input  temperature  profiles  (dotted  line)  compared  with  the  true  profiles 
and,  as  before,  the  computed  heat  fluxes. 

Although  the  method  is  relatively  insensitive  to  the  specific  value  of  3, 
care  should  be  taken  to  insure  that  it  is  of  the  right  order  of  magnitude.  To 
illustrate  this,  the  first  case  for  low  noise  was  repeated  with  3  much  to  large 
(1)  and  S  much  too  small  (10-®).  These  results  are  shown  in  Figures  5a  and  5b 
respectively. 


16 


NSWC  TR  82-32 


5.  CONCLUDING  REMARKS 

The  method  described  in  this  report  provides  an  effective  tool  for  computing 
heat  fluxes  from  thermocouple  data  in  wind  tunnel  experiments.  For  three  typical 
cases  with  characteristic  experimental  noise  levels  excellent  agreement  with  the 
known  solution  was  obtained.  The  method  is  not  completely  automated  and  requires 
user  judgement  in  its  application.  The  user  must  specify  a  S  value  based  on  the 
expected  noise  level  and  often  a  trial  and  error  procedure  is  required  to  optimize 
results.  Also  the  regularization  technique  leads  to  degradation  of  results  as 
t  -♦  T.  This  problem  can  be  circumvented  by  recording  data  past  the  en>J  of  the  run. 
Finally,  in  the  straightforward  extension  to  the  nonlinear  case  and  multidimensional 
problems  the  method  becomes  extremely  cumbersome  and  costly.  However,  the  general 
approach  seems  promising  and  further  work  may  make  it  possible  to  circumvent  the 
above  difficulties. 
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APPENDIX  A 

COMPUTER  PROGRAM  LISTING 

program  heatinv ( input , output , tape 5= input ,  tape6=output , 
tapell,tapel2) 

this  program  computes  heat  fluxes  on  the  front  face 
from  temperature  data  measured  at  the  back  face, 
analysis  of  several  sets  of  data  is  possible  in 
each  run. 

the  data  from  the  back  face  is  read  off  of  tapell 
each  set  of  data  is  assumed  to  be  written  unformatted 
with  the  number  of  points  first  and  then  the  data 
values,  it  is  assumed  that  there  are  the  same  number 
of  points  in  each  data  set.  otherwise,  the  matrix 
must  be  reformed. 

the  program  outputs  the  computed  values  of  qdot  (the 
heat  fluxes)  are  output  on  tapel2  for  post  processing, 
see  subroutine  output  for  a  description  of  the  write 
statements 

input  of  relevant  parameters  is  specified  in  a  namelist, 
dimension  ip(100) 

common  a (100 , 100) , alpha ,alam,dx ,dt , ulf (100) ,g(100) ,b(100) 

, f ( 101 ) ,jj,nt,uinit,aa (3,401) ,u(401) , rhs (401) ,beta , tmax ,fac 

data  stored  in  blank  common: 

a  stores  the  matrix  which  approximates  the 

action  of  the  regularized  operator  which 
computes  backface  temperatures  from  given 
fluxes  on  the  front  face 

aa  stores  the  matrix  used  in  solving  the  heat 

equation 

g  stores  the  vector  which  specifies  the  matrix 

ulf  stores  the  measured  values  on  the  left  face 

f  stores  the  fluxes  that  are  to  be  used  when 

solving  the  heat  equation 

u  stores  the  temperatures  profiles  while  solving  the 

heat  equation 

rhs  right  hand  side  for  solving  the  heat  equation 
dx  value  of  delta  x  used  to  solve  the  heat  equation 

dt  value  of  delta  t  used  to  solve  the  heat  equation 

alam  parameter  used  for  crank  nicholson 

uinit  initial  value  of  the  temperature 
tmax  the  length  of  the  run  in  seconds 

fac  ratio  of  u  sub  x  to  qdot 
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c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 

c 
c 
c 
c 
c 
c 
c 

j j  =40 1  $  uinit=0.  $  ipr  =  l  $  beta=l.e-4 

read(5, input) 
wri te (6 , input) 
mdim=100 

dx=l ./float ( j j-1) 
dt=l./float(nt) 

alam=alpha*144 . *tmax/delta**2*dt/2./dx**2 
fac=12.*ak/delta 
call  setmat 

call  decomp(nt,mdim,a , ip) 
do  100  i=l,nthc 
read ( 11) ndp, (ulf (n) , n=l , ndp) 
if (ndp.ne .nt)  go  to  1000 
i f (ndp.gt .mdim) go  to  1001 
call  setrhs 

call  solve (nt ,mdim,a fb , ip) 
do  300  n=l,nt 
f (n+1) =b (ip(n) ) 

300  continue 
f  (l)-0. 

call  output { ipr , i ) 

100  continue 
stop 

1000  write (6 , 1002) i 

1002  format(*  number  of  data  points  for  thermocouple  *fi5, 
1  *  is  inconsistent  with  the  input  value*) 

stop 

1001  wr i te ( S , 1003) mdim , i 

1003  format(*  there  are  more  than  *,i5,*  data  points  for  * 

1  , ‘thermocouple  *,i5) 


nt  number  of  data  points  and  number  of  steps  taken 

in  solving  the  heat  equation 

beta  this  is  the  regularization  parameter.  beta  should 
be  selected  so  that  the  estimated  12  norm 
printed  in  output  correspond  to  the  estimated 
rms  error  in  the  measured  data,  when  the  12  norm 
is  smaller  than  the  estimate  results  can  become 
oscillatory,  (remedy:  increase  beta)  if  larger 
the  results  can  be  degraded,  best  value  of 
beta  is  when  these  two  rms  values  agree, 
the  default  value  of  beta  is  suitable  for  the 
expected  noise  levels  in  wind  tunnel  data  reduction 

namel is t/ input/ j  j , ipr ,nt , nthc ,beta,uinit,alpha,ak, tmax , delta 

data  specified  by  the  namelist 

j j ,nt , beta ,uini t , alpha ,ak , del ta  are  specified  above 

alpha  thermal  diffusivity  in  ft**2  per  sec 

ak  thermal  conductivity  in  btu  per  (ft  sec  deg  f) 

delta  thickness  of  the  plate  in  inches 
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end 

subroutine  setmat 


NSWC  TR  82-32 


c  this  routine  computes  the  matrix  corresponding  to 
c  the  normal  equations  o£  the  regularized  approximation 
c  to  the  integral  equation 

common  a (100, 100) , alpha ,alam,dx ,dt,ulf (100) ,g(100) ,b(100) 

1  #£(101) ,jj,nt,uinit,aa(3,401) ,u(401) ,rhs(401) ,beta , tmax , f ac 

do  100  n=l,nt 
f (n)*0. 

100  continue 

f  (2) =1.  $  f (nt+l) =0. 

call  heat(g) 
do  120  id  *  1,  nt 
n  =»  nt  -  id  +  1 
do  110  jd  *  1,  id 
m  *  nt  -  jd  +  1 
a (m,n)  =  g( jd) *g(id) 

if(jd  .gt.  1)  a(m,n)  =  a(m,n)  +  a(m+l,n+l) 
a(n,m)  =  a(m,n) 

110  continue 
120  continue 

here  we  add  the  stabilization  part  of  the  matrix 

scrl  *  2.*beta*(dt/3.+l./dt) 
scr2  *  beta*(dt/6.  -  l./dt) 
ntml  =  nt  -  1 
do  200  n  =  1,  ntml 
a(n,n)  =  a(n,n)  +  scrl 
a(n+l,n)  =  a(nrn+l)  =»  scr2+a  (n,n+l) 

200  continue 

a(nt,nt)  *  a (nt ,nt) +scr 1/2 . 
return 
end 

subroutine  setrhs 

c  this  routine  computes  the  right  hand  side 
c  of  the  normal  equations  from  the  right  hand 
c  side  of  the  least  squares  problem 
c 

common  a(100,100) , alpha, alam,dx,dt,ulf (100) ,g(100) ,b(100) 

1  ,f(101) ,jj,nt,uinit,aa(3,401) ,u(401) ,rhs(401) , beta , tmax , f ac 

do  100  i=l,nt 
b ( i ) *0  . 
do  120  j»i,nt 
b(i)=b(i)+ulf ( j) *g( j-i+1) 

120  continue 
100  continue 
return 
end 
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subroutine  heat(culf) 

this  routine  computes  the  solution  of  the 
c  heat  equation  using  crank  nicholson 

c  the  fluxes  are  specified  at  x=l  in  the 
c  the  vector  f  and  the  bacface 

c  temperatures  are  returned  in  the  vector  culf 

c 

dimension  culf(l) 

common  a ( 100 , 100) , alpha , al am ,dx ,dt,ulf(100) ,g ( 100 ) ,b (100 ) 

1  , f (101) ,jj,nt,uinit,aa(3,401) ,u(401) , rhs (401) , beta , tmax , f ac 

do  10  j*l,jj 
u ( j ) *uini t 
10  continue 
call  form 
j jml=j j-1 
do  100  n*l,nt 

rhs (1)  = (1.-2. *alam) *u ( 1) +2 . *alam*u ( 2) 

rhs(jj)=(l.-2.*alam) *u ( j j ) +2.*alam* (  u ( j jml ) +dx* ( f (n) +f (n+1) ) ) 
do  110  j=2,jjml 

rhs ( j ) * (l.-2.*alam) *u( j ) +alam* (u ( j-1 ) +u ( j+1 ) ) 

110  continue 

call  tr islv (aa , rhs , j j ) 
do  120  j*l,jj 
u( j)=rhs( j) 

120  continue 

culf (n) =u ( 1 ) 

100  continue 
return 
end 

subroutine  form 
c 

c  this  routine  forms  the  matrix  used  in  advancing 
c  the  heat  equation 

c 

common  a (100 , 100) , alpha ,alamfdx,dt,ulf(100) ,g(100) ,b(100) 

1  , f ( 101) ,jj,nt,uinit,aa(3,401) ,u(401) ,rhs(401) ,beta, tmax , fac 

do  100  j=l,jj 
aa(2,j)al.+2.*alam 
aa ( 1 , j ) =aa  ( 3 , j) =-alam 
100  continue 

aa (1, j j)*aa ( 3 , 1) =-2 . *alam 
call  tr idem (aa , j j ) 
return 
end 

subroutine  tridcm(a,jj) 
c 

c  decomposes  a  tridiagonal  matrix  stores  in  band  form 
c 

dimension  a(3,l) 
do  100  j*2,jj 
jml* j-1 

a (1, j)*scr*a (1, j)/a (2, jml) 
a  (2, j) =a  (2,  j) -scr*a(3, jml) 

100  continue 
return 
end 
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subroutine  trislv (a,b, j j) 
c 

c  uses  result  of  tridcm  to  solve  a  x  =  b 

c  storing  the  result  in  b 

c 

dimension  a(3,l),b(l) 
do  100  j=2,jj 
b( j)=b(j)-a(l, j)*b( j-1) 

100  continue 

b(jj)=b(jj)/a(2,jj) 
do  200  jd=2,jj 
j~ j  j+2- jd 
jml= j-1 

b  ( jml  )=*(b(jml)-a(3,jml)  *b( j) )/a  (2,  jml) 

200  continue 
return 
end 

subroutine  output(ipr , i) 
c 

c  computes  output  for  each  thermocouple 

c 

common  a(100,100) , alpha ,alam,dx ,dt , ul f (100) ,g(100 ) ,b(100) 

1  , f ( 101 ) ,jj,nt,uinit,aa (3,401) ,u  (401 ) , rhs (401) ,beta , tmax ,  f ac 

write (12) nt , (f loat (n) *dt*tmax , f (n+1 ) *f ac ,n=l ,nt) 
if ( ipr .eq. 0) return 
write(6,1003) i 

1003  format (*lanalysis  for  thermocouple  *,i5) 
call  heat(b) 
write (6,1000) 

1000  format(*Otime*,llx,*backface  temperatures ( f) *, 6x , 

1  ,15x,*  heat  flux  (btu/ft* , 2h** ,*  sec)*/15x,*  computed* ,7x, 

2  *measured*,7x,*differences*//) 
sum=0 . 

do  100  n=l,nt 
t* float ( n) *tmax*dt 
sum=sum+ (b (n) -ulf (n) ) **2 

write (6, 1001) t,b (n) ,ulf (n) ,b (n) -ulf (n) ,f (n+1) *fac 

1001  format(lx,lp5el5.7) 

100  continue 

sum=sqrt (sum/ float (nt) ) 
write  (6,1002)sum 

1002  format(*  the  12  norm  of  differences  in  the  back  face* 

1  ,*  temperatures  is  *,lpel5.7) 

return 

end 
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